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1  Abstract 


In  this  work,  we  carry  out  optimized  spread-spectrum  concealment  of  multiuser 
data  under  a  given  digital  image  (or  audio,  or  video  sequence).  First,  the  overall 
image/host  medium  is  pre-processed  into  transform-domain  small  blocks  from 
which  host  vectors  are  obtained  via  zig-zag  scanning  vectorization.  Multiuser 
data  hiding  is  performed  in  the  generated  host  vectors.  Under  this  data  hiding 
system  model,  we  calculate  an  orthogonal  set  of  embedding  spread-spectrum 
signatures  that  achieves  maximum  sum  signal-to-interference-plus-noise  ratio 
(sum-SINR)  at  the  output  of  the  linear-filter  receivers  for  any  fixed  embedding 
amplitude  values.  Then,  for  any  given  total  embedding  distortion  constraint, 
we  present  the  optimal  multi-signature  assignment  and  amplitude  allocation 
that  maximizes  the  sum  capacity  of  the  concealment  procedure.  The  prac¬ 
tical  implication  of  the  reported  results  is  sum-SINR,  sum-capacity  optimal 
multiuser/multi-signature  spread-spectrum  data  hiding  in  the  digital  medium. 
Extensive  experiments  that  we  carried  out  demonstrate  the  effectiveness  of  the 
proposed  methods. 

We  also  consider  the  problem  of  extracting  blindly  data  embedded  over  a  wide 
band  in  a  spectrum  (transform)  domain  of  a  digital  medium  (image,  audio, 
video).  We  develop  a  novel  multi-carrier/signature  iterative  generalized  least- 
squares  (M-IGLS)  core  procedure  to  seek  unknown  data  hidden  in  hosts  via 
multi-carrier  spread-spectrum  embedding.  Neither  the  original  host  nor  the  em¬ 
bedding  carriers  are  assumed  available.  Experimental  studies  on  images  show 
that  the  developed  algorithm  can  achieve  recovery  probability  of  error  close  to 
what  may  be  attained  with  known  embedding  carriers  and  host  autocorrelation 
matrix. 
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2  Introduction 


Embedding  digital  data  in  digital  host  images  has  a  broad  range  of  applications, 
from  annotation  and  single-stream  media  merging  (text /audio/image)  to  water¬ 
marking  and  covert  communications  [1] ,  [2]  when  embedding  is  to  be  carried  out 
in  a  manner  imperceptible  by  the  human  eye. 

There  are  different  approaches  in  the  literature  to  embed  digital  data  in 
a  host  image.  Simple  direct  hiding  in  a  given  host  is  considered  in  [3]- [6]. 
Embedding  information  after  full  frame  discrete  fourier  transform  (DET)  or 
discrete  cosine  transform  (DCT)  image  transformation  is  discussed  in  [7]-[10]. 
Concealing  data  after  block  DET  or  DCT  transformation  of  the  host  image 
is  presented  in  [11], [12].  Concealing  data  after  wavelet  transformation  of  the 
original  host  image  is  described  in  [13], [14]. 

In  this  work,  we  consider  the  emerging  concept  of  multiuser  (multi-signature) 
data  embedding  where  multiple  messages  are  hidden  with  different  embedding 
spread-spectrum  signatures  in  the  host  image  [15], [16].  The  theoretical  chal¬ 
lenges  of  multiuser  spread-spectrum  data  embedding  in  part  parallel  those  prob¬ 
lems  encountered  in  the  field  of  code-division  multiple-access  (CDMA)  com¬ 
munications  [17].  Eor  CDMA  signature  design,  in  the  theoretical  context  of 
complex/real- valued  signature  sets,  the  early  work  of  Welch  [18]  on  total-square- 
correlation  (TSC)  bounds  was  followed  up  by  direct  minimum-TSC  designs 
[19], [20]  and  iterative  distributed  optimization  algorithms  [21], [22].  New  bounds 
on  the  TSC  of  binary  signature  sets  were  found  [23]  that  led  to  minimum-TSC 
optimal  binary  signature  set  designs  for  almost  all  signature  lengths  and  set 
sizes  [23]-[25].  The  sum  capacity,  total  asymptotic  efficiency,  and  maximum 
squared  correlation  of  the  minimum-TSC  binary  sets  were  evaluated  in  [26]. 
New  bounds  and  optimal  designs  for  minimum  TSC  quaternary  signature  sets 
are  derived  in  [27].  The  problems  of  periodic  and  aperiodic  TSC  are  treated  in 
[28], [29].  Binary /quaternary  adaptive  signature  design  over  multipath  channels 
is  discussed  in  [30], [31]. 

In  this  present  work,  based  on  the  optimal  orthogonal  carriers  that  max¬ 
imize  the  sum  capacity  of  the  multiple-access  colored  Gaussian  vector  chan¬ 
nel  presented  in  [32],  we  give  the  orthonormal  set  of  signatures  for  multiuser 
spread-spectrum  data  hiding  in  digital  images  that  offers  maximum  sum  signal- 
to-interference-plus-noise  ratio  (sum-SINR)  embedding  in  arbitrary  transform 
domain  images  with  any  given  embedding  amplitude  values.  Moreover,  for  any 
given  total  host  distortion  budget  we  present  a  power  (amplitude)  allocation 
scheme  that  maximizes  the  Shannon  sum-capacity  of  the  multiuser  data  embed¬ 
ding  system  under  the  assumption  of  Gaussian  distributed  (transform-domain) 
host  data.  These  theoretical  findings  establish  optimality  of  the  Gkizeli-Pados- 
Medley  multisignature  eigen-design  algorithm  [16]  under  the  general  require¬ 
ment  of  an  orthogonal  multiuser  signature  set.  Experimental  studies  and  com¬ 
parisons  included  herein  illustrate  the  theoretical  results. 

The  notation  used  in  this  report  is  as  follows:  j-}^  denotes  the  transpose 
operation;  M"  denotes  the  n  dimensional  real  field;  E{-}  represents  statistical 
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expectation;  I„  is  the  identity  matrix  of  size  n  x  n.  We  use  boldface  lower¬ 
case  letters  to  denote  column  vectors  and  boldface  uppercase  letters  to  denote 
matrices. 


3  System  Model 

Consider  a  host  image  H  €  ,  where  ^  is  the  image  pixel  alphabet  and 

X  N2  is  the  image  size  in  pixels.  The  overall  host  image  H  is  partitioned 
into  P  small  blocks.  Each  small  block  has  N1N2/P  pixels.  Under  the  multiuser 
data  embedding  model  introduced  in  [16],  each  small  block  Hi,  H2, . . . ,  Hp  is 
about  to  transport  K  hidden  digital  information  bits,  one  for  each  different  user 
potentially.  We  will  perform  data  hiding  in  the  real  two-dimensional  transform 
domain 

After  transform  calculation  for  each  small  block,  the  matrices  of  the  transform- 
domain  coefficients  p=  1,  2, . . . ,  P,  are  vectorized  by  conventional  “zig¬ 

zag  scanning”  to 

Vec{.^(Hp)}  =  [l7(Hp)i,i  l7(Hp)2.i  ^(Hp)i,2 

.^(Hp)3.l  ^(Hp)2.2  ^(Hp)i,3  ...]^  (1) 

where  ^{'H.p)ij  denotes  the  {i,j)th  element  of  matrix  .T(Hp).  Note  that 
Vec{.^(Hp)}  e  p  =  1, 2,  •  •  •  ,  P. 

The  final  host  vectors 


XpGR^,  p=l,2,...,P,  (2) 

of  length  L  <  N1N2/P  are  formed  directly  from  any  subset  of  coefficients  of 
Vec{l3^(Hp)},  p  =  1,  2,  •  •  •  ,  P.  For  example,  we  can  have  host  vector  length  L  = 
N1N2/P  —  1  by  excluding  only  the  dc  coefficient  since  modification 

of  i^(Hp)i^i  is  known  to  lead,  in  general,  to  visible  image  change  in  the  pixel 
domain.  We  show  the  diagram  of  generating  the  host  vectors  from  a  given  host 
image  in  Fig.  1. 

The  autocorrelation  matrix  of  the  transform-domain  host  vectors  x  is  defined 
as 

Rx=E  {xx^}  = 

It  is  easy  to  verify  that  for  general  natural  images  R^,  ^  uIl  ,  a  >  0;  that  is,  Ra, 
is  not  constant-valued  diagonal  or  “white”  in  field  language.  The  host  image 
example  of  Elaine  is  shown  in  Fig.  2  where  =  {0, 1,  •  •  •  ,  255}^^®^^^®. 

The  image  is  partitioned  into  P  =  =  1024  8  x  8  small  blocks.  After 

two-dimensional  DOT  transformation  of  the  small  blocks  and  zig-zag  scanning 
vectorization,  we  obtain  Vec{.i?^(Hi)},  Vec{.i5^(H2)},  •  •  • ,  Vec{.T(Hio24)}  with 
individual  vector  length  64.  We  exclude  only  the  dc  coefficients  to  form  the 
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final  host  vectors  xi,  X2,  •  •  • ,  X1024  of  length  L  =  63.  The  host  data  autocorre¬ 
lation  matrix  R^,  is  shown  in  Fig.  3.  We  can  see  that  the  host  vectors  can  be 
considered/act  as  colored  noise  to  data  to  be  embedded. 

Multiuser  data  hiding  can  carried  out  in  the  host  vectors  xi,  X2,  •  •  • ,  xp. 
Assuming  that  there  are  K  users  of  interest  or  signatures  for  spread-spectrum 
data  embedding,  the  host  vectors  are  modified  to 

K 

y  =  ^  AibiSi  -f  X  +  n  (4) 

i=l 


where  bi,b2,  ■  ■  ■  ,bK  G  {±1}  are  the  individual  message  bits  embedded  simulta¬ 
neously  in  X  €  {xi,  X2,  •  •  • ,  Xp}  with  corresponding  amplitudes  Ai  >  0  and 
normalized  signatures  ||  ||=  1,  i  =  1,  2, . . . ,  K.  The  additive  vector 

n  ^  pT(0,CT^Ip)  accounts  in  the  model  for  possible  external  white  Gaussian 
noise  of  variance  Embedded  bits,  host  vectors  x,  and  noise  vectors  n  are 
considered  to  be  statistically  independent  from  each  other.  The  embedded  bits 
themselves  are  considered  as  Bernoulli  probability- 1/2  random  variables  that 
are  independent  across  time  and  users/messages. 

For  independent  embedded  bits  (or  orthogonal  signatures) ,  we  can  calculate 
the  mean-square  distortion  over  the  original  host  image  as 


^  =  E 


K 


(5) 


The  distortion  caused  by  each  individual  embedded  message  z,  i  = 
is 


9,  =  E\\\A,hs£]=Al 


(6) 


Assume  now  that  given  y,  message  j  €  {1,  2, . . . ,  AT}  is  the  message  of  inter¬ 
est.  With  signal  of  interest  AjbjSj  and  respective  total  disturbance 
X  -|-  n  from  (4) ,  the  linear  filter  that  operates  on  y  and  offers  maximum  SINK 
at  its  output  can  be  calculated  using  the  Cauchy-Schwarz  inequality  [33]  to  be 


'WmaxSiNRj  =  arg  max 


E 


( 


:  AibiSi  -I-  X  -f  n 


=  R,.^s,- 
h  ^ 

where  matrix  R/^  is  defined  as 


R/.  =  E 


K 


s,-  H-  X  +  n 


AibiS^  +  X  -f  n 


Ra;  -I-  (T^Il  +  y]  AjsiS^ 


(7) 


(8) 
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Note  that  R/j  is  a  function  of  Si,  •  •  •  ,  Sj_i,  Sj+i,  •  •  •  ,sk  and  independent  of  Sj. 

Then,  the  output  SINK  value  when  we  use  the  filter  'WmaxSiNRj  can  be 
calculated  as 

SINRj  =  A^sJ  I  Ra;  +  ctAl  +  ^  Sisf  j  sj  (9) 

\  / 

=  A^jsjRjhj.  (10) 

We  define  the  metric  of  sum-SINR  for  the  multiuser  data  embedding  system  as 

K 

sumSINR  =  SINRj.  (11) 

1=1 

As  mentioned  in  [16],  eq.  (36),  we  can  independently  maximize  SINRi 
for  user  1  with  respect  to  Si;  then,  maximize  SINR2  for  user  2  with  respect 
to  S2,  etc.  After  one  optimization  cycle  over  Si,  S2,  •••  ,S/f,  we  can  update 
R/i,  R/2,  •  •  ■  I'^/K  accordingly  and  run  a  second  cycle  over  Si, S2,  •  •  •  , S/f .  We 
may  continue  with  optimization  and  update  cycles  until  numerical  convergence  is 
observed.  Regretfully,  there  is  no  known/guaranteed  optimality  of  this  described 
scheme.  In  this  work,  we  will  present  a  one-shot  sum-SINR  (and  sum-capacity) 
optimal  signature  set  assignment  that  operates  directly  on  the  transform-domain 
host  data  autocorrelation  matrix  R^,  of  (3). 


4  Optimal  Multisignature  Embedding 


First,  we  recall  the  Matrix  Inversion  Lemma  [34]  (also  known  as  Woodbury’s 
Identity) , 

(B  +  UCV)-i  =  -  B”iU(C-i  -I-  VB~iU)-^VB-\  (12) 


where  B,  U,  C  and  V  all  denote  matrices  of  appropriate  size. 

Then,  we  begin  a  tedious  -yet  all  important-  algebraic  derivation  of  the 
SINR  expression  for  user  j,  j  =  1,2,-  -  •  ,Ar,  in  (9).  The  application  of  the 
Matrix  Inversion  Lemma  on  equation  (9)  leads  to 


SINRj 


K 


A|sJ  Rj;  -I-  cr^Ii  -1-  Y,  A 


•4N  ( 

(Rx  +  ctAl  +  ^ 

( 

(r®  + 

-1 


SINRf~^'>  -  tY~^\ 


j<K-l, 
j  =  K 

(13) 
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where 


SINd^~^^  =  J 


-1 


Z,j  , 


j<K-l, 

A'jsJ  (^Ra-  +  CT^Ii  +  Sj, 

j  =  K, 


(14) 


and 


rp(K-l)  A 
3 


(K-l)  A 

Pi  =  •{ 


_ _ 

j<K-l, 

_ _ 

iM?f_i+s]^_i(Rx+<T2ii+E(il"  Afsisf  ’ 

j  =  K, 


[  sj  (Ro;  +  Cr^Ii  +  A%sJ )  , 

j<K-l, 


(15) 


(16) 


sj  (R-o;  +Cr2lL  Sk_i, 

j  =  K. 

Using  again  the  Matrix  Inversion  Lemma  on  SINR^^~^^  in  (14),  we  obtain 
analogous  expressions  SINRj^~^\  ■  Continuing  on  repetitively, 

we  obtain  SINR^^~^\  etc.,  of  the  general  form  SINr!'J^\ 

and  l<n<it'  —  1,  given  below 


SINR^r^  =  I 


1 .4’sj  (r, +aUi. +2r=u* /ifsisr) 


( 


3  <  n, 

Rx  +  ct^Il  +  YJlZi  A^Sisf )  Sj , 


(17) 


3  >n, 


rpin)  A  . 

j  ~  ' 


Adp<")|^ 


lM^  +  l+sJ+l(R-x+<T2li+5:;"^j_.^j.  A^^Sisf)  S„  +  1  ■ 

3  <  n, 

_ _ 

i/Al+sT{R^+anL+T.?A^  ’ 

3  >  n, 


(18) 


sj  +  Er=l,i/J  S„+1, 

j  <  n, 


(n)  A 

P  =  ^ 


T 


(^Ra;  +  (J^Il  +  I]r=/  SisJ 


-1 

Sni 

J  >  n. 


(19) 


By  (9),  after  K  —  1  application  of  the  Matrix  Inversion  Lemma  we  reach 
SINK,  = 

J  J 


K-l 

=  SINRf^  -  rj”\  (20) 

n—1 

where 

SINRf'^  =  (R,  +  s,.  (21) 

Accordingly,  we  can  calculate  the  sum-SINR  metric  as 

K 

sumSINR  =  SINRj 

i=i 

K  K  K-l 

tY-  (22) 

j=l  j=l  n=l 

Let  {qi,  q2, . . . ,  Ql}  be  the  L  eigenvectors  of  Ra,  with  corresponding  eigen¬ 
values  Ai  >  A2  >  . . .  >  Al-  The  proof  of  the  following  Lemma  comes  directly 
from  [32]. 


Lemma  1  With  orthonormal  signature  sets  {si,  S2, . . . ,  K  <  L,  for  mul¬ 
tiuser  spread- spectrum  data  embedding  and  corresponding  fixed  embedding  ampli¬ 
tudes  Al  >  A2  >  ■  ■  ■  >  Ak  >  0,  SINR^^^  is  maximized  to  Yi)+(t'^ 

when  Si,  S2,  sk  are  assigned  as  the  K  smallest- eigenvalue  eigenvectors  of 
the  image  host  vectors  autocorrelation  matrix  Ra;,  i.e.,  Sj  =  j  = 

1,2, . . . ,  K .  At  the  same  time,  when  Sj  =  qi_(j_i),  j  =  1,2, . . . ,  K ,  =  0 

for  every  j  =  1,2, ...  ,K  and  n  =  1, . . . ,  K  —  1.  □ 


Equipped  with  the  result  of  Lemma  1  on  equation  (22),  we  return  to  the 
problem  of  optimal  multiuser  data  embedding  in  host  images  by  (4)  and  con¬ 
sider  the  multiuser  performance  metrics  sumSIN R  =  'Yj=i  SINRj  and  sum 
capacity  Csurm  defined  as  the  maximum  information  that  can  be  obtained  about 
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the  hidden  variables  bi, . . .  ,bK  by  examining  the  received  vector  y  over  all  pos¬ 
sible  input  probability  distributions  /(&i, .  • . ,  6_fc), 

Csum=  max  l{bi,...,bK-,y)-  (23) 

If  the  transform-domain  host  data  x  can  be  assumed  as  Gaussian,  the  sum 
capacity  of  the  data  embedding  scheme  is  [16] 

Csun.  =  ^  log2  det  [1l  +  (R.  +  a^W-^SA^S'^]  (24) 


where  S  =  [si,  S2, . . . ,  is  the  L  x  K  matrix  formed  by  the  spread-spectrum 

embedding  signatures  and  A  =  diag(Ai,  A2, . . . ,  Ak)  is  the  diagonal  matrix  of 
the  embedding  amplitudes. 

The  following  theorem,  built  on  Lemma  1,  establishes  the  optimal  orthonor¬ 
mal  multiuser  data  embedding  signature  set  for  host  images  under  fixed  ampli¬ 
tude  embedding. 


Theorem  1  (Optimal  Multi-signature  Embedding) 

For  orthonormal  signature  sets  {si,  S2, . . . ,  },  K  <  L,  for  multiuser  spread- 

spectrum  data  embedding  and  corresponding  fixed  embedding  amplitudes  Ai  > 
4I2  >  .  • .  >  Ak  >  0,  the  sum-SINR  is  maximized  to  sumSINR-max  = 

J2j=i  \)+cr'^  when  Si,  S2,  . . .,  sk  are  assigned  as  the  K  smallest- eigenvalue 
eigenvectors  of  the  image  host  vectors  autocorrelation  matrix  Ra,,  i.e.,  Sj  = 


Ql — (j — 1 )  ?  J  1 7  2 , ... ,  A . 

If  the  transform- domain  host  data  x  are  Gaussian  distributed,  the  same  signa¬ 
ture  assignment  maximizes  sum  capacity  to  {Csum)max  ~ 


1 

2  =  l 


log2  ( 


1 


bits  per  K-symbol  embedding. 


□ 


Hence,  by  Theorem  1,  we  proved  that  the  sum-capacity  optimality  condition 
presented  as  necessary  for  multiuser  spread-spectrum  data  hiding  in  Lemma  1 
of  [16]  is  also  sufficient  for  general  orthonormal  signature  set  design.  We  can 
further  consider  individual  embedding  amplitudes  as  design  parameters  as  well, 
within  a  total  distortion  budget.  In  other  words,  we  can  search  for  the  optimal 
amplitude  assignment  that  maximizes  sum  capacity  for  Gaussian  host  vectors 
subject  to  a  total  allowable  distortion  constraint  ^  The  jointly 

optimal  power  allocation  and  signature  assignment  is  presented  in  Theorem  2 
below. 


Theorem  2  (Optimal  Multisignature  Embedding  and  Power  Alloca¬ 
tion) 

For  orthonormal  signature  sets  and  a  given  total  embedding  distortion  budget 
the  optimal  ( signature,  amplitude )  pairs  that  maximize  the  sum  capacity  for 
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multiuser  data  embedding  in  (transform- domain)  Gaussian  hosts  are 
(sj  =  qL-(j-i),  =  [-  {Xl-H-i)  +cr^)  +m]^)  > 

where  [a:]’*'  =  max{x,0)  and  /i  is  the  Kuhn-Tueker  [35]  coefficient  chosen  such 
that  the  distortion  constraint  '3  =  equality.  □ 

The  amplitude/power  allocation  method  described  in  Theorem  2  can  be 
viewed  as  an  eigen  domain  “waterfilling”  procedure  [35]. 


5  Experimental  Studies  on  Embedding 

In  this  section,  we  will  first  experiment  on  one  image  as  an  example  and  then 
show  average  results  over  the  entire  USC-SIPI  database  [36]  of  44  miscellaneous 
images. 

We  first  use  as  a  host  the  familiar  grayscale  “Boat”  image  of  size  512x512 
in  Fig.  4(a).  We  partition  the  image  into  small  blocks  of  size  8x8  and  calculate 
the  DCT  transform  of  each  block.  Then,  we  remove  the  dc  coefficient  from 
the  zig-zag  scanned  DCT-domain  vectors  and  create,  this  way,  a  host  set  of 
512^/8^  =  4096  vectors  of  length  L  =  8x8— 1  =  63.  We  hide  K  =  15  data 
messages  via  multiuser  spread-spectrum  embedding  (15  bits  in  each  block,  each 
bit  on  a  different  signature)  and  include  also  additive  white  Gaussian  noise  of 
variance  cr^  =  3dB.  Fig.  4(b)  shows  the  Boat  image  after  signature-set  optimal 
multiuser  spread-spectrum  data  hiding  by  Theorem  1  of  this  report  at  1^  =  20dB 
total  distortion.  The  individual  message  amplitudes/distortions  are  fixed  at 
3i,  3i  =  3\  —  IdB,  •  •  • ,  ^15  =  ^14  —  IdB  {IdB  decrease  for  each  successive 
message).  Fig.  4(c)  shows  the  Boat  image  after  multiuser  spread-spectrum 
data  hiding  at  =  20dB  total  distortion  with  jointly  optimal  signature  set 
and  amplitudes  by  Theorem  2  of  this  report.  No  perceptual  difference  can  be 
observed  by  naked  eye  between  the  three  images  in  Fig.  4. 

In  Fig.  5,  we  plot  the  sum-SINR  of  the  messages  reception  process  when 
the  total  embedding  distortion  3  varies  from  12  to  32dB  and  embedding  is 
carried  out  with  either  (i)  arbitrary  or  (ii)  optimal  signatures  by  Theorem  1. 
For  this  example,  the  gain  in  sum  SINK  by  the  use  of  the  optimal  signature 
set  of  Theorem  1  is  at  least  5dB  and  grows  as  the  total  allowed  distortion 
increases.  To  translate  sum-SINR  to  the  more  immediate  and  familiar  metric 
of  average  bit-error-rate,  in  Fig.  6  we  plot  the  corresponding  average  bit  error 
rate  (BER)  of  the  recovered  concealed  messages.  We  can  observe,  for  example, 
that  when  we  use  the  optimal  set  of  Theorem  1,  to  recover  hidden  messages 
with  bit  error  rate  10“^  we  need  to  cause  only  about  31dB  distortion  (including 
the  accounted  white  noise).  To  have  similar  recovery  error  rate  with  arbitrary 
sequence-set  embedding  is  not  even  practically  possible  per  Fig.  6. 

In  Fig.  7,  we  examine  the  sum  capacity  of  multiuser  spread-spectrum  data 
hiding  under  (i)  arbitrary  signature  set  design,  (ii)  signature  optimization  alone 
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by  Theorem  1,  and  (iii)  optimal  signature  and  power  allocation  by  Theorem  2. 
At  32dB  total  distortion,  an  optimized  signature  set  (with  fixed  per  message 
distortion  at  IdB  increments)  offers  a  gain  of  about  10  bits  per  embedding  at¬ 
tempt  over  arbitrary  signature  sets.  About  3  more  bits  are  added  when  signature 
optimization  is  combined  with  optimal  power  allocation. 

As  a  concluding  evaluation  effort,  we  carry  out  the  experiments  of  Figs.  5 
and  7  over  the  entire  USC-SIPI  database  [36]  of  44  miscellaneous  images  shown 
in  Fig.  8  (16  color  and  28  monochrome;  14  of  size  256  x  256,  26  of  size  512  x  512, 
and  4  of  size  1024x  1024.)  Fig.  8  shows  average  sum  SINK  versus  total  distortion 
over  the  whole  database.  Fig.  9  shows  average  sum  capacity  results  over  the 
database.  The  average  database  findings  of  Figs.  9  and  10  are  quite  similar 
to  the  individual  Boat  results  in  Figs.  5  and  7  offering  credible  experimental 
support  for  the  optimized  multiuser  embedding  procedures  described  in  this 
report. 


6  Conclusions  on  Embedding 

We  considered  the  problem  of  multiuser  spread-spectrum  data  hiding  in  transform- 
domain  hosts  (images  in  particular,  herein)  and  identified  the  orthonormal  sig¬ 
nature  set  that  offers  maximum  sum  SINK  embedding  for  any  fixed  embed¬ 
ding  amplitude  values.  We  showed  that  the  set  is  also  sum  capacity  optimal  in 
terms  of  bits  per  multiuser  embedding  under  the  assumption  that  the  transform- 
domain  host  data  is  Gaussian.  When  there  is  flexibility  in  assigning  amplitudes 
across  users  under  a  total  host  distortion  constraint,  we  derived  the  user  am¬ 
plitude  values  that  not  only  meet  the  total  constraint  but  also  maximize  sum 
capacity. 
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Figure  1:  Host  vectors  generation  diagram. 
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Figure  2:  Host  image  example  of  Elaine. 
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Figure  3:  Host-vector  autocorrelation  matrix  of  Elaine. 
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(c) 


Figure  4:  (a)  Original  Boat  image  example  (512x512  grayscale),  (b)  Boat 
image  after  optimal  multi-signature  embedding  {K  =  15  messages  of  size  4096 
bits  each,  total  distortion  20dB,  fixed  per  message  distortion  =  ^i  —  IdB, 
i  =  1,  ■  ■  ■  ,  K  —  1,  =  2>dB).  (c)  Boat  image  after  optimal  signature  and 

power  allocation  {K  =  15  messages  of  size  4096  bits  each,  total  distortion  20dB, 
al  =  MB). 
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Figure  5:  Sum  SINK  versus  total  distortion  (Boat  image,  K 
9i  -  UB,  i  =  I,--  -  ,i4:  -  1,  0-2  =  MB). 


BER 


Figure  6:  BER  versus  total  distortion  (Boat  image,  K  =  15,  ^i+i  =  ^i  —  IdB, 
i  =  I,--  -  ,R:  -  1,  cr2  =  MB). 
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Figure  8:  USC-SIPI  image  database  of  44  images. 


23 


Total  distortion  D  (dB) 


Figure  9:  Sum  SINK  versus  total  distortion  (average  findings  over  USC-SIPI 
image  database  [36],  K  =  15,  ^i+i  =  —  IdB,  i  =  1,  -  ■  ■  ,  K  —  1,  =  MB). 
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Figure  10:  Sum  capacity  versus  total  distortion  (average  findings  over  USC-SIPI 
image  database  [36],  K  =  15,  =  idB). 


7  The  Extraction  Problem 


Digital  data  embedding  in  digital  media  is  an  information  technology  field  of 
rapidly  growing  commercial  as  well  as  national  security  interest.  Applications 
may  vary  from  annotation,  copyright-marking,  and  watermarking,  to  single¬ 
stream  media  merging  (text,  audio,  image)  and  covert  communication  [1].  In 
annotation,  secondary  data  are  embedded  into  digital  multimedia  to  provide  a 
way  to  deliver  side  information  for  various  purposes;  copyright-marking  may 
act  as  permanent  “iron  branding”  to  show  ownership;  fragile  watermarking 
may  be  intended  to  detect  future  tampering;  hidden  low-probability-to-detect 
(LPD)  watermarking  may  serve  as  identification  for  confidential  data  validation 
or  digital  fingerprinting  for  tracing  purposes  [2]-[4].  Covert  communication  or 
steganography,  which  literally  means  “covered  writing”  in  Greek,  is  the  process 
of  hiding  data  under  a  cover  medium  (also  referred  to  as  host),  such  as  im¬ 
age,  video,  or  audio,  to  establish  secret  communication  between  trusting  parties 
and  conceal  the  existence  of  embedded  data  [5]- [9].  As  a  general  encompassing 
comment,  different  applications  of  information  hiding,  such  as  the  ones  iden¬ 
tified  above,  require  different  satisfactory  tradeoffs  between  the  following  four 
basic  attributes  of  data  hiding  [10]:  (i)  Payload  -  information  delivery  rate;  (ii) 
robustness  -  hidden  data  resistance  to  noise/disturbance;  (Hi)  transparency  - 
low  host  distortion  for  concealment  purposes;  and  (iv)  security  -  inability  by 
unauthorized  users  to  detect /access  the  communication  channel. 

Recently,  developing  data  embedding  technologies  are  being  seen  to  pose  a 
threat  to  personal  privacy,  commercial,  and  national  security  interests  [11],  [12]. 
In  this  work,  we  focus  our  attention  on  the  blind  recovery  of  secret  data  hidden 
in  medium  hosts  via  multi-carrier/signature  direct-sequence  spread-spectrum 
(DS-SS)  transform  domain  embedding  [13]-[20].  Neither  the  original  host  nor 
the  embedding  carriers  (signatures  or  spreading  sequences)  are  assumed  known 
(fully  blind  data  extraction).  This  blind  hidden  data  extraction  problem  has 
also  been  referred  to  as  “Watermarked  content  Only  Attack”  (WOA)  in  the 
watermarking  security  context  [21]-[24]. 

While  passive  detection-only  of  the  presence  of  embedded  data  is  being  inten¬ 
sively  investigated  in  the  past  few  years  [25]- [33],  active  hidden  data  extraction 
is  a  relatively  new  branch  of  research.  In  blind  extraction  of  SS  embedded  data, 
the  unknown  host  acts  as  a  source  of  interference/disturbance  to  the  data  to 
be  recovered  and,  in  a  way,  the  problem  parallels  blind  signal  separation  (BSS) 
applications  as  they  arise  in  the  fields  of  array  processing,  biomedical  signal 
processing,  and  code-division  multiple-access  (CDMA)  communication  systems 
[34]- [38].  Under  the  assumption  that  the  embedded  secret  messages  are  inde¬ 
pendent  identically  distributed  (i.i.d.)  random  sequences  and  independent  to 
the  cover  host,  independent  component  analysis  (ICA)  may  be  utilized  to  pur¬ 
sue  hidden  data  extraction  [24],  [39].  However,  ICA-based  BSS  algorithms  are 
not  effective  in  the  presence  of  correlated  signal  interference  as  is  the  case  in 
SS  multimedia  embedding  and  degrade  rapidly  as  the  dimension  of  the  carrier 
(signature)  decreases  relative  to  the  message  size.  In  [19],  an  iterative  general¬ 
ized  least  squares  (IGLS)  procedure  was  developed  to  blindly  recover  unknown 
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messages  hidden  in  image  hosts  via  SS  embedding.  The  algorithm  has  low  com¬ 
plexity  and  strong  recovery  performance.  However,  the  scheme  is  designed  solely 
for  single- carrier  SS  embedding  where  messages  are  hidden  with  one  signature 
only  and  is  not  generalizable  to  the  multi- carrier  case.  Realistically,  an  em- 
bedder  would  favor  multi-carrier  SS  transform-domain  embedding  to  increase 
security  and/or  payload  rate. 

In  this  work,  we  develop  a  novel  multi-carrier  iterative  generalized  least 
squares  (M-IGLS)  algorithm  for  SS  hidden  data  extraction  that,  to  the  best  of 
the  authors’  knowledge,  appears  for  the  first  time  in  the  broad  communication 
theory  and  systems  literature.  For  improved  recovery  performance,  in  particular 
for  small  hidden  messages  that  pose  the  greatest  challenge,  experimental  studies 
indicate  that  a  few  independent  M-IGLS  re-initializations  and  executions  on 
the  host  can  lead  to  hidden  data  recovery  with  probability  of  error  close  to 
what  may  be  attained  with  known  embedding  carriers  and  known  original  host 
autocorrelation  matrix.  Applications  of  the  developed  algorithm  are,  of  course, 
not  limited  to  attacking  steganographic  covert  communications  by  recovering 
the  secret  embedded  messages.  Since  the  carriers  are  also  jointly  estimated 
with  the  embedded  data,  the  developed  scheme  can  also  be  used  for  complete 
message  removal  or  tampering  attack  as  well  by  reinserting  a  fabricated  message 
in  place  of  the  original.  From  the  opposite  data  embedding  point  of  view,  the 
developed  algorithm  can  be  treated  as  a  tool  to  test  security  robustness  of  SS 
data  hiding  schemes. 

The  following  notation  is  used  below.  Boldface  lower-case  letters  indicate 
column  vectors  and  boldface  upper-case  letters  indicate  matrices;  M  denotes  the 
set  of  all  real  numbers;  denotes  matrix  transpose;  Tr{-}  is  matrix  trace;  II 
is  the  L  X  L  identity  matrix;  sgn{-}  denotes  zero-threshold  quantization;  and 
E{-}  represents  statistical  expectation.  Finally,  /  I,  ||  •  ||,  and  |j  •  jj^.  are  the  scalar 
magnitude,  vector  norm,  and  matrix  Frobenius  norm,  respectively. 

8  Multi-carrier  SS  Embedding  and  Extraction: 
Problem  Formulation 

Gonsider  a  host  image  H  €  -v^here  Ai  is  the  finite  image  alphabet  and 

Ni  X  N2  is  the  image  size  in  pixels.  Without  loss  of  generality,  the  image  H 
is  partitioned  into  M  local  non-overlapping  blocks  of  size  ■  Each  block. 
Hi,  H2, ....,  Hm,  is  to  carry  K  hidden  information  bits  {KM  bits  total  image 
payload).  Embedding  is  performed  in  a  2-D  transform  domain  T  (such  as  the 
discrete  cosine  transform,  a  wavelet  transform,  etc.).  After  transform  calculation 
and  vectorization  (for  example  by  conventional  zig-zag  scanning),  we  obtain 
T(H^)  G  M  M  ,  m  =  1,  2, . . . ,  M.  From  the  transform  domain  vectors  T{lim) 
we  choose  a  fixed  subset  of  L  <  coefficients  (bins)  to  form  the  final  host 
vectors  x(m)  €  m  =  1,2, ,  M.  It  is  common  and  appropriate  to  avoid 
the  dc  coefficient  (if  applicable)  due  to  high  perceptual  sensitivity  in  changes  of 
the  dc  value. 

The  autocorrelation  matrix  of  the  host  data  x  is  an  important  statistical 
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quantity  for  our  developments  and  is  defined  as  Rx  =  E{xx^}  =  X]m=i  x(m)x(m)^. 
It  is  easy  to  verify  that  in  general  Rx  7^  Q;Il,  a  >  0;  that  is,  Rx  is  not  constant- 
value  diagonal  or  “white”  in  field  language.  For  example,  8x8  DCT  with  63-bin 
host  data  formation  (excluding  only  the  dc  coefficient)  for  the  256  x  256  gray¬ 
scale  Baboon  image  in  Fig.  11(a)  gives  the  host  autocorrelation  matrix  Rx  in 
Fig.  11(b)  [20]. 

8.1  Multi-carrier  SS  Embedding 

We  consider  K  distinct  message  bit  sequences,  {6^(1),  6^(2), ...,  5fc(M)},  k  = 

1,2, . . . ,  K,  bkim)  e  {±1},  m  =  each  of  length  M  bits.  The  K 

message  sequences  may  be  to  be  delivered  to  K  distinct  corresponding  recipients 
or  they  are  just  K  portions  of  one  large  message  sequence  to  be  transmitted 
to  one  recipient.  In  particular,  the  mth  bit  from  each  of  the  K  sequences, 
bi (m) , . . .  ,bK {m) ,  is  simultaneously  hidden  in  the  mth  transform-domain  host 
vector  x(m)  via  additive  SS  embedding  by  means  of  K  spreading  sequences 
(carriers)  ||sfc||  =  1,  /p  =  1,  2, . . . ,  K, 

K 

y(m)  =  ^  Akhk{m)sk  +  x(m)  +  n(m),  m  =  1,  2, . . . ,  M,  (25) 

k^l 

with  corresponding  amplitudes  Ak  >  Q,  k  =  1, . . .  ,K .  For  the  sake  of  generality, 
n(m)  represents  potential  external  white  Gaussian  noise^  of  mean  0  and  autocor¬ 
relation  matrix  >  0.  It  is  assumed  that  bk{m)  behave  as  equi-probable 

binary  random  variables  that  are  independent  in  m  (message  bit  sequence)  and 
k  (across  messages).  The  contribution  of  each  individual  embedded  message  bit 
bk  to  the  composite  signal  is  A}.bk^k  and  the  block  mean-squared  distortion  to 
the  original  host  data  x  due  to  the  embedded  k  message  alone  is 

Vk=E{\\AkSkbkf}  =  Al,  k=l,2,...,K.  (26) 

Under  statistical  independence  of  messages,  the  block  mean-squared  distortion 
of  the  original  image  due  to  the  total,  multi-message,  insertion  of  data  is  T>  = 

The  intended  recipient  of  the  fcth  message  with  knowledge  of  the  fcth  carrier 
Sfe  can  perform  embedded  bit  recovery  by  looking  at  the  sign  of  the  output  of 
the  minimum-mean-square-error  (MMSE)  hlter  'WMMSE,k  =  Ry 

bk{m)  =  sgni^liMSE  ,fcy(w)}  =  sgn{s^R- V(w)}  (27) 

where  Ry  is  the  autocorrelation  matrix  of  the  host-plus-data-plus-noise  vectors 

K 

—  IE{yy^}  =  Rx  +  ^  AlskS^  +  (28) 

k^l 

^Additive  white  Gaussian  noise  is  frequently  viewed  as  a  suitable  (most  entropic)  model 
for  general  quantization  errors,  channel  transmission  disturbances,  and/or  image  processing 
attacks  [40]. 
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The  autocorrelation  matrix  Ry  can  be  estimated  by  sample  averaging  over  the 
set  of  M  received  vectors  {y(TO)}^^;^,  Ry  =  Using  Ry 

in  (27)  in  place  of  Ry,  we  obtain  what  is  known  as  the  sample-matrix-inversion 
MMSE  (SMI-MMSE)  detector  implementation  [34], 


8.2  Formulation  of  the  Extraction  Problem 

To  blindly  extract  spread-spectrum  embedded  data  from  a  given  host  image, 
the  analyst  needs  first  to  convert  the  host  to  observation  vectors  of  the  form 
of  y(m),  m  =  1, . . .  ,M,  in  (25).  This  requires  knowledge  of  (i)  the  partition, 
(n)  transform  domain,  (in)  subset  of  coefficients,  and  (iv)  number  of  carriers 
used  by  the  embedder.  The  host  image  partition  (and  block  size  N1N2/M  in 
our  notation)  may  be  estimated  by  neighboring-pixels  difference  techniques  as 
in  [30].  Regarding  the  subset  of  coefficients  used  in  embedding,  the  conserva¬ 
tive  approach  is  to  assume  that  all  coefficients  are  used,  except  maybe  the  dc 
value,  and  set  accordingly  L  =  N1N2/M  —  1.  The  number  of  carriers  K  can 
be  estimated  by  SS  signal  population  identification  algorithms  such  as  in  [41]. 
Finally,  determination  of  the  transform  domain  used  in  embedding  seems  to  be 
a  hurdle  not  yet  tackled  by  current  research.  The  natural  approach  would  be 
to  consider  individually  and  exhaustively  one  transform  at  a  time  starting  from 
the  most  common  (for  example,  2D-DCT,  common  wavelet  transforms,  and  so 
on). 

In  this  report,  we  focus  the  technical  presentation  solely  after  the  point  that 
the  analyst  obtains  transform-domain  observations  in  the  form  of  y(m)  in  (25), 
upon  performing  appropriate  image  partition  and  transform  calculation.  We 
denote  the  combined  “disturbance”  to  the  hidden  data  (host  plus  noise)  by 
z(m)  =  x(to)  -|-  n(m)  and  rewrite  SS  embedding  by  (25)  as 

K 

y(m)  =  ^Akbk{ni)sk  +  z(m),  m  =  1, . . .  (29) 

k^l 

where  z(m)  is  modeled  as  a  sequence  of  zero- mean  (without  loss  of  generality) 
vectors  with  autocovariance  matrix  R^  =  E{zz^}  =  Rx-l-cr^I.  Let  =  A^Sk  G 
k  =  1, . . . ,  K,he  the  amplitude-including  embedding  carriers.  Then,  we  can 
further  reformulate  SS  embedding  as 

K 

y(m)  =  +  z(m)  (30) 

k^l 

=  Vb(m) -I- z(to),  m  =  1, . . . ,  M,  (31) 

where  V  =  [vi, . . . ,  v^^]  G  is  the  amplitude- including  carrier  matrix  and 

b(TO)  G  {±1}^^^  is  the  vector  of  bits  embedded  in  the  mth  host  block.  For 
notational  simplicity,  we  can  write  the  whole  observation  data  in  the  form  of 
one  matrix 

Y  =  VB  -b  Z  (32) 
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where  Y  =  [y(l)  y(2)  ...  y(M)]  G  B  =  [b(l)  b(2)  ...  b(M)]  G  {±1}^^’^ , 

and  Z  =  [z(l)  z(2)  ...  z(M)]  G 

Our  objective  is  to  blindly  extract  the  unknown  hidden  data  B  from  the 
observation  matrix  Y  without  prior  knowledge  of  the  embedding  carriers 
and  amplitudes  Ak,  /c  =  1, . . . ,  if,  in  V  =  [AiSi, . . . ,  Ak&k]  or  the  host  medium 
itself  x(l), . . . ,  x(M)  in  Z  =  [x(l)  +  n(l), . . . ,  x(M)  +  n(M)]. 


9  Hidden  Data  Extraction 


If  Z  were  to  be  modeled  as  Gaussian  distributed,  the  joint  maximum-likelihood 
(ML)  estimator  of  V  and  decoder  of  B  would  be 

V,B  =  arg  min  ||R;^  (Y  -  VB)||2  (33) 

VeR-^  X  K 


where  multiplication  by  Rz  ^  can  be  interpreted  as  prewhitening  of  the  com¬ 
pound  observation  data.  If  Gaussianity  of  Z  is  not  to  be  invoked,  then  (33)  can 
be  simply  referred  to  as  the  joint  generalized  least-squares  (GLS)  solution^  of 

V  and  B. 

The  global  GLS-optimal  message  matrix  B  in  (33)  can  be  computed  inde¬ 
pendently  of  V  by  exhaustive  search  over  all  possible  choices  under  the  criterion 
function  |jRz  ^YRbIIf, 

B  =  arg  min  IjR^^YRBlI?.  (34) 

where  Rb  =  I  —  B^(BB^)“^B.  The  derivation  of  (34)  is  provided  in  the 
Appendix.  Exhaustive  search  has,  of  course,  complexity  exponential  in  KM 
(total  size  of  hidden  messages  in  bits).  We  consider  this  cost  unacceptable  and 
attempt  to  reach  a  quality  approximation  of  the  solution  of  (34)  (or  (33),  to 
that  respect)  by  alternating  generalized  least-squares  estimates  of  V  and  B, 
iteratively,  as  described  below. 

Pretend  B  is  known.  The  generalized  least-squares  estimate  of  V  is 
VoLs  =  arg  min  ||Rz”"  (Y  -  VB)||2 

veKixK 

=  YB^(BB'^)-A  (35) 


Pretend,  in  turn,  that  V  is  known.  Then,  the  least-squares  estimate  of  B  over 
the  real  field  is 


B 


real 


arg  min  ||Rz  ^  (Y  -  VB)||2 

(V^R-iV)-iV^R-iY. 


(36) 


^Generalized  least-squares  solutions  are  weighted  least-squares  (WLS)  solutions  with  op- 

_  I 

timal  weighting  matrices,  here  Rz  ^  ,  that  yield  the  lowest  variance  of  the  estimation  error 
[35],  [36]. 
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Observing  that 


(V^R-iV)-iV^R-i  =  (V^R-iV)-iV'^R-\  (37) 

we  rewrite 

=  (V^R;1V)-1V^R;1Y  (38) 

and  suggest  the  approximate  binary  message  solution 

=  arg  min  HR,”* (Y  -  VB)||2 
cs  sgn{(V^R;iV)-iV^R;iY}.  (39) 

The  proofs  of  (35),  (36),  and  (37)  are  provided  in  the  Appendix. 

The  multi-carrier  iterative  generalized  least-squares  (M-IGLS)  procedure  sug¬ 
gested  by  the  two  equations  (35)  and  (39)  is  now  straightforward.  Initialize  B 
arbitrarily  and  alternate  iteratively  between  (35)  and  (39)  to  obtain  at  each 
step  conditionally  generalized  least  squares  estimates  of  one  matrix  parameter 
given  the  other.  Stop  when  convergence  is  observed.  Notice  that  (39)  utilizes 
knowledge  of  the  autocorrelation  matrix  Ry,  which  can  be  estimated  by  sample 
averaging  over  the  received  data  observations,  Ry  =  ^  y(^)y(''^)^-  The 

M-IGLS  extraction  algorithm  is  summarized  in  Table  1.  Superscripts  denote  it¬ 
eration  index.  The  computational  complexity  of  each  iteration  of  the  M-IGLS 
algorithm  is  0{2K^  -\-  2LMK  -\-  K^{3L  -\-  M)  -\-  L'^K)  and,  experimentally,  the 
number  of  iterations  executed  is  between  20  and  50  in  general. 

For  the  sake  of  mathematical  accuracy,  we  recall  that  in  least-squares  there 
is  always  a  symbol  sign  (phase  in  complex  domains)  ambiguity  when  joint  data 
extraction  and  carrier  estimation  is  pursued  (i.e.,  data  bits  G  {±1}^  on 
carrier  G  have  the  same  least-squares  error  with  data  bits  — on  carrier 
— Sfc,  k  =  1, . . .  ,K).  The  sign-ambiguity  problem  can  be  overcome  with  a  few 
known  or  guessed  data  symbols  for  supervised  sign  correction^  [42].  Moreover, 
in  a  multi-carrier  least-squares  scenario  as  the  one  that  we  face  herein,  the  index 
association  remains  unresolved  (i.e.,  given  a  recovered  (message,  carrier)  pair 
(b,  s),  the  corresponding  index  A:G{l,...,Ar}in  (25)  cannot  be  obtained).  To 
the  extend  that  the  application  of  the  work  presented  in  this  report  is  to  simply 
extract  blindly  the  embedded  bits  with  the  least  possible  errors,  the  carrier 
indexing  problem  is  not  dealt  with  any  further. 


Returning  to  the  proposed  data  extraction  algorithm,  we  understand  that 
with  arbitrary  initialization  convergence  of  the  M-IGLS  procedure  described  in 
Table  1  to  the  optimal  GLS  solution  of  (33)  is  not  guaranteed  in  general.  Exten¬ 
sive  experimentation  with  the  algorithm  in  Table  I  indicates  that,  for  sufficiently 
long  messages  hidden  by  each  carrier  (M  =  4Kbits  or  more,  for  example),  sat¬ 
isfactory  quality  message  decisions  B  can  be  directly  obtained.  However,  when 

®If  the  embedded  data  are  encrypted,  then  both  options  and  — bj.  must  be  separately 
decrypted  and  investigated  for  sign  correction  for  each  message  k  =  I, K. 
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Table  1:  Multi-carrier  Iterative  generalized  least-squares  Data  Extraction 


1) 

d:= 

0;  initialize  B*^°) 

G  {±l}^fxM 

arbitrarily. 

2) 

d:= 

d  -f  1; 

yid) 

:=  Y(B(‘^-i))^ 

(BG-i))(bG 

-..)T] 

1 

5 

bG) 

:=sgn  j((V('^)) 

^  (VW)^RyiY|. 

3) 

Repeat  Step  2  until  B^'^)  = 

the  message  size  is  small,  M-IGLS  may  very  well  converge  to  inappropriate 
points/solutions.  The  quality  (generalized-least-squares  fit)  of  the  end  conver¬ 
gence  point  depends  heavily  on  the  initialization  point  and  arbitrary  initializa¬ 
tion  -which  at  first  sight  is  unavoidable  for  blind  data  extraction-  offers  little  as¬ 
surance  that  the  iterative  scheme  will  lead  us  to  appropriate,  “reliable”  (close  to 
minimal  generalized  least-squares  fit)  solutions.  To  that  respect,  re-initialization 
and  re-execution  of  the  M-IGLS  procedure,  say  P  times,  is  always  possible.  To 
assess  which  of  the  P  returned  solutions,  say  (Vi,  Bi), . . . ,  (Vp,  Bp),  has  supe¬ 
rior  generalized-least-squares  fit,  we  simply  feed  (V^B^)  to  (33)  (using  Ry  in 
place  of  Rz)  and  choose 

Vflnai,Bfi„ai  =  arg  min  HR^"  (Y  -  VB)||^  (40) 

(V3)e{(Vi,Bi),....(Vj,,Bp)} 

The  computational  complexity  of  the  P-times  re- initialized  M-IGLS  is,  of  course, 
0(PD{2K^  +  2LMK  +  K^{3L  +  M)  +  Li^K))  where  D  represents  the  number 
of  internal  iterations  in  d  in  Table  1. 


10  Experimental  Studies  on  Extraction 

A  technically  firm  and  keen  measure  of  quality  of  a  hidden-message  extraction 
solution  is  the  difference  in  bit-error-rate  (BER)  experienced  by  the  intended 
recipient  and  the  analyst.  The  intended  recipient  in  our  studies  may  be  using  any 
of  the  following  three  message  recovery  methods:  (i)  Standard  carrier  matched- 
filtering  (ME)  with  the  known  carriers  s^,  k  =  1,...,A';  (ii)  sample-matrix- 
inversion  MMSE  (SMI-MMSE)  filtering  with  known  carriers  Sj,  and  estimated 
host  autocorrelation  matrix  Ry  (see  (27));  and  (in)  ideal  MMSE  filtering  with 
known  carriers  Sj,  and  known  true  host  autocorrelation  matrix  Rx,  which  serves 
as  the  ultimate  performance  bound  reference  for  all  methods.  In  terms  of  blind 
extraction  (neither  nor  Rx  known),  we  will  examine:  (iv)  The  developed 
M-IGLS  algorithm  in  Table  1  with  P  =  20  re-initializations  and,  for  comparison 
purposes,  the  performance  of  two  typical  independent  component  analysis  (IGA) 
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based  blind  signal  separation  (BSS)  algorithms  (v)  FastICA  [44],  and  (vi)  JADE 
[45]. 

We  first  consider  as  a  host  example  the  gray-scale  512  x  512  “Baboon”  image. 
We  perform  8x8  block  DCT  embedding  by  (25)  over  all  bins  except  the  dc 
coefficient  with  K  =  A  distinct  arbitrary  carriers  k  =  1,. . .  ,K.  The 

hidden  message  embedded  by  each  carrier  is  =  4,  096  bits  long.  The  per- 
message  block  mean  square  distortion  due  to  each  embedded  message  is  set  to  be 
the  same  for  all  messages,  i.e.  Vk  =  A\  =  fc  =  1, . . . ,  4.  With  per- message  8 x 
8-block  MSE  distortion  Vk,  k  =  1, . . . ,  AT,  the  peak  signal-to-noise  ratio  (PSNR) 
of  the  image  due  to  embedding  can  be  calculated  by  PSNR  A  201ogig(255)  — 
101ogig(^^;^  I?fe/64).  Another  metric  that  reflects  the  relationship  between 
host  and  embedding  distortion  is  the  block  document-to-watermark  power  ratio 
(DWR)  defined  as  DWR  A  lOlog^Qd^  —  101ogjQ(^^^  1?^)  where  cr^  A  Tr{Rx} 
is  the  (total)  host  block  variance.  The  value  of  depends  on  the  nature  of  each 
host  image  and  is  provided  in  each  experiment  that  we  run  (see  figure  captions) 
to  facilitate  translation  by  the  reader  between  MSE  and  DWR  if  desired.  For 
the  sake  of  generality,  in  our  studies  we  also  incorporate  white  Gaussian  noise 
of  variance  =  3dB. 

Fig.  12  shows  the  average  BER  (over  all  AT  =  4  messages)  of  all  methods 
(i)  through  (vi)  listed  above  as  a  function  of  the  host  block  distortion  per- 
message.  FastICA  and  JADE  have  computational  complexity  0{2(K  —  1){K  + 
M)  +  5MK{K  +  l)/2)  per  iteration  and  0(K(K  -  1)(AK^  +  21K  +  75)/2), 
respectively.  In  particular,  on  an  Intel  i5-2550K  3.40GHz  processor  running 
standard  Matlab  software  for  experimentation,  the  average  execution  time  of 
the  M-IGLS  algorithm  with  P  =  20  initializations  was  1.51  sec,  the  average 
execution  time  of  FastICA  was  0.20  sec,  and  the  average  execution  time  of 
JADE  was  0.08  sec.  While  the  two  independent /principal-component  methods 
(FastICA  and  JADE)  are  failing  to  carry  out  effective  hidden  data  extraction,  to 
our  satisfaction  M-IGLS  analysis  is  very  close  in  BER  performance  to  the  ideal 
MMSE  detector  bound  in  which  both  the  embedding  carriers  and  the  clean  host 
autocorrelation  matrix  Rx  are  treated  as  perfectly  known. 

In  Fig.  13,  we  repeat  the  exact  same  experimental  study  on  the  smaller 
256  X  256  version  of  the  Baboon  image  in  Fig.  11(a)  with  AT  =  4  hidden 
messages  of  length  only  ^||-  =  1,  024  bits  per  message  (compared  to  4, 196 
bits  per  message  in  Fig.  12).  Comparing  with  Fig.  12,  the  gap  between  M- 
IGLS  and  ideal  MMSE  increases  as  the  hidden  message  size  decreases,  but 
the  extraction  performance  of  M-IGLS  can  still  be  deemed  satisfactory.  For 
additional  experimental  validation,  the  studies  of  Fig.  12  and  Fig.  13  are 
repeated  on  the  familiar  “Boat”  image  (shown  in  Fig.  14)  in  its  512  x  512  and 
256  X  256  gray-scale  versions  (Fig.  15  and  Fig.  16,  correspondingly). 

To  examine  the  behavior  of  M-IGLS  under  increasing-density  small-message 
hiding,  we  consider  the  256  x  256  gray-scale  “F-I6  Aircraft”  image  (shown  in 
Fig.  17)  with  AT  =  4  and  K  =  8  hidden  messages  of  length  1Kbit  each.  The 
recovery  performance  plots  for  A"  =  4  and  AT  =  8  are  given  in  Figs.  18  and  19, 
correspondingly. 
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An  encompassing  conclusion  over  all  executed  experiments  is  that  M-IGLS 
remains  a  most  effective  technique  to  blindly  extract  hidden  messages,  while 
extraction  becomes  more  challenging  as  the  length  of  the  hidden  message  per 
used  embedding  carrier  decreases  or  the  number  of  hidden  messages  (number 
of  used  carriers)  increases.  It  is  also  worth  pointing  out  that,  in  these  ex¬ 
perimental  studies,  M-IGLS  may  outperform  (in  moderate  to  high  distortion 
values)  SMI-MMSE  in  which  the  true  carriers/signatures  are  known.  This  is 
because  SMI-MMSE  suffers  from  performance  degradation  due  to  small-sample- 
support  adaptation  (estimation  of  matrix  Ry).  The  unsatisfactory  performance 
of  the  ICA-based  methods  is  due  to  the  interference  from  high- amplitude  (low- 
frequency)  host  coefficients.  To  demonstrate  this  point,  in  Fig.  20  we  repeat 
the  exact  same  experiment  of  Fig.  12  using  this  time  only  the  L  =  20  highest- 
frequency  DCT  coefficients  as  our  host  vector.  It  can  be  observed  that,  in  this 
moderate  host  interference  environment,  ICA-based  methods  can  provide  satis¬ 
factory  performance  (not  superior  to  M-IGLS,  however).  Of  course,  we  may  not 
expect  that  data  are  always  embedded  exclusively  in  low-amplitude  coefficients 
alone. 

Next,  for  the  sake  of  enhanced  experimental  credibility,  we  examine  the 
average  performance  of  the  proposed  M-IGLS  algorithm  over  a  large  image 
database.  The  experimental  image  data  set,  [46]  and  [47]  combined,  consists 
of  more  than  11,500  8-bit  gray-scale  photographic  images  which  have  great 
variety  (e.g.,  outdoor/indoor,  daylight /night,  natural/man-made)  and  different 
sizes.  We  embed  one  up  to  five  messages,  K  G  {1,2,...,  5},  via  multi-carrier 
SS  embedding  with  arbitrary  carriers  and  payload  between  0.016  and  0.078 
bits  per  pixel  (bpp).  The  length  of  the  embedding  carriers  varies  between  30 
and  63,  L  G  {30, 31, . . . ,  63}.  Recovery  performance  plots  are  given  in  Fig. 
21.  Similar  conclusions  can  be  drawn  as  in  the  previous  individual  image  host 
experimentations. 

While  our  blind  data  extraction  algorithmic  development  was  based  on  the 
most  common  SS  embedding  form  (1)  for  convenience  in  presentation,  the  de¬ 
veloped  algorithm  can  also  be  applied  to  more  advanced  SS  embedding  schemes 
such  as  improved  spread-spectrum  (ISS)  [13]  and  correlation-aware  improved 
spread-spectrum  (CAISS)  [43].  In  Fig.  22,  we  go  again  over  the  whole  [46],  [47] 
databases  under  ISS  embedding  and  in  Fig.  23  under  GAISS  embedding  (with 
amplitude-proportion  parameter  rj  =  0.7)^.  It  can  be  noted  from  Figs.  22,  23 
that  M-IGLS  analysis  can  also  carry  out  effective  hidden  data  extraction  for  the 
ISS  and  GAISS  schemes. 


11  Conclusions  on  Extraction 

We  considered  the  problem  of  blindly  extracting  unknown  messages  hidden  in 
image  hosts  via  multi-carrier/signature  spread-spectrum  embedding.  Neither 

■^Both  ISS  [13]  and  CAISS  [43]  are  proposed  as  single  carrier  embedding  schemes  {K  =  1 
in  the  experiments). 
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the  original  host  nor  the  embedding  carriers  are  assumed  available.  We  devel¬ 
oped  a  low  complexity  multi-carrier  iterative  generalized  least-squares  (M-IGLS) 
core  algorithm.  Experimental  studies  showed  that  M-IGLS  can  achieve  prob¬ 
ability  of  error  rather  close  to  what  may  be  attained  with  known  embedding 
signatures  and  known  original  host  autocorrelation  matrix  and  presents  itself  as 
an  effective  countermeasure  to  conventional  SS  data  embedding/hiding^. 


®In  [39],  Bas  and  Cayre  presented  an  interesting  signature-based  additive  embedding  ap¬ 
proach  different  to  (25)  that  is  host-vector-by-host-vector  dependent  and  would  withstand 
IGLS-based  analysis.  The  embedding  is,  however,  sensitive  to  noise  which  would  lead  to 
high  recovery  error  rates  by  intended  recipients  and  limit  the  applicability  to  general  covert 
communication  problems. 
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APPENDIX 

A.  Derivation  of  (34) 

The  minimization  in  (33)  can  be  carried  out  in  two  steps.  First,  we  minimize 
(33)  with  respect  to  V:  V  =  YB^(BB^)“^  (see  also  Appendix  B).  Then, 
substituting  V  back  into  (33)  we  obtain 

B  =  arg  min  llRz  ®(Y- YB^(BB'^)-1B)||2  (41) 

=  arg  min  ||Rz  ®  Y(I  -  B^(BB'^)-1B)||2  (42) 

Be{±i}^x“ 

=  arg  min  HR^^YRbII'  (43) 

where  Rb  =  I-B^(BB'^)-iB.  ■ 


B.  Proof  of  (35) 

The  GLS  cost  function  in  (33)  can  be  rewritten  as 

J  ^  ||Rz”^Y-R;5vB|1|,  (44) 

=  Tr  |r-1YY'^|  -  Tr  |r- AB^V^|  -  Tr  |r- ABY^|  +  Tr  |r- ABB^ylb^ 

where  Tr{-}  denotes  the  trace  of  a  matrix. 

For  a  given  message  matrix  B,  the  GLS  optimal  estimate  of  V  can  be  obtain 
by  differentiating  the  cost  function  J  with  respect  to  and  setting  the  outcome 
equal  to  the  zero  matrix, 

^  =  -R^-^YB^  +  R7iV(BB^)  =  0  ^  V  =  YB^(BB^)-i.  (46) 


C.  Proof  of  (36) 

We  rewrite  the  GLS  cost  function  in  (45)  as 

J  =  Tr  |r-1YY'^|  -  Tr  jv^R"  AB^|  -  Tr  {R-^VBY^}  +  Tr  {V^R^^VBB^}  .  (47) 

Pretend  that  V  is  known  and  relax  the  domain  of  the  symbol  information  matrix 
to  the  real  space,  B  G  The  GLS  optimal  estimate  of  B  G  can 

be  calculated  again  by  differentiation 

^  =  -V'^R-^Y  +  V^R^-^VB  =  0  ^  B  =  (V^R-^V)~^V'^R-^Y.  (48) 
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D.  Proof  of  (37) 

Since  Ry  =  Elyy"^}  =  VV"^  +  R„  by  the  Matrix  Inverse  Lemma  (also  known 
as  Woodbury’s  matrix  identity)  [48],  we  obtain 

R;1  =  R-i  -  R-iV(I  +  V^R-iV)-iV^R-i.  (49) 

Then, 

V^R-^V  =  V^R-^V- V^R-iV(I  +  V^R7iV)-iV^R-iV 
=  v^R-iv[i  -  (I  +  v'^R-^vy^v^R-^v] 

=  V^R-iV(I  + V^R-iV)-i[(I  + V'^R-^V)  -  V^R-^V] 

=  V^R-iV(I  + V^R-^V)-!.  (50) 

By  the  property  of  the  inverse  of  a  product  of  matrices  [48],  the  inverse  of 
(V^R-^V)  is 

(V^R-^V)-!  =  (I  +  V^R7iV)(V^R-^V)-i 

=  (V^R-iV)-i+I.  (51) 

We  combine  the  results  of  (49)  and  (51)  and  finally  obtain 

(V^R”iV)-iV^R-i  =  ((V'^R-iV)-i  +  I)  (R7I  -  R-iV(I  +  V^R^Vy^V'^Ry) 

=  (V^R-iV)-iV^R-\  (52) 
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Figure  11:  (a)  Baboon  image  example  H  €  {0, 1, (b)  Host  data 
autocorrelation  matrix  (8x8  DCT,  63-bin  host)  [20]. 
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Figure  12:  Average  BER  versus  per-message  block  distortion  (512  x  512  Baboon, 
L  =  63,  AT  =  4  messages  of  4Kbits  each,  cr^  =  3dB,  =  46.49dB). 


Figure  13:  Average  BER  versus  per-message  block  distortion  (256  x  256  Baboon, 
L  =  63,  AT  =  4  messages  of  1Kbit  each,  =  3dB,  cr^  =  45.45dB). 


(a)  (b)  (c) 

Figure  14:  512  x  512  gray-scale  Boat  image. 


Figure  15:  Average  BER  versus  per-message  block  distortion  (512  x  512  Boat, 
L  =  63,  AT  =  4  messages  of  4Kbits  each,  cr^  =  3dB,  cr^  =  44.15dB). 
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Figure  16:  Average  BER  versus  per-message  block  distortion,  256  x  256  Boat, 
L  =  63,  AT  =  4  messages  of  1Kbit  each,  =  3dB,  cr^  =  45.57dB). 


Figure  17:  256  x  256  gray-scale  Aircraft  image. 
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Figure  18:  Average  BER  versus  per-message  block  distortion  (256  x  256  F16 
Aircraft,  L  =  63,  K  =  4:  messages  of  1Kbit  each,  =  3dB,  =  46.23dB). 


Figure  19:  Average  BER  versus  per-message  block  distortion  (256  x  256  E16 
Aircraft,  L  =  63,  K  =  8  messages  of  1Kbit  each,  cr^  =  3dB,  =  46.23dB). 


Figure  20:  Average  BER  versus  per-message  block  distortion  (512  x  512  Baboon, 
L  =  20  highest-frequency  coefficients,  K  =  A  messages  of  4Kbits  each,  = 
3dB,  al  =  46.49dB). 


Figure  21:  Average  BER  versus  per-message  block  distortion  (average  findings 
over  a  dataset  of  more  than  11,  500  images  [46],  [47],  K  =  1,  L  G  {30, 31, ... ,  63}, 
cr^  =  3dB,  average  cr^  =  41.63dB). 


Figure  22:  Average  BER  versus  per-message  block  distortion  (average  findings 
over  a  dataset  of  more  than  11,  500  images  [46],  [47],  ISS  embedding  [13],  AT  =  1, 
L  G  {30,  31, ... ,  63},  cr^  =  3dB,  average  =  41.63). 


Figure  23:  Average  BER  versus  per-message  block  distortion  (average  findings 
over  a  dataset  of  more  than  11,500  images  [46],  [47],  CAISS  embedding  [43], 
(r;  =  0.7),  a:  =  1,  L  e  {30, 31, ... ,  63},  al  =  3dB,  average  cr^  =  41.63dB). 
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